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Abstract 

We discuss efficient solutions to systems of shifted linear systems 
arising in computations for oscillatory hydraulic tomography (OHT). 
The reconstruction of hydrogeological parameters such as hydraulic 
conductivity and specific storage, using limited discrete measurements 
of pressure (head) obtained from sequential oscillatory pumping tests, 
leads to a nonlinear inverse problem. We tackle this using the quasi- 
linear geostatistical approach [13J. This method requires repeated so- 
lution of the forward (and adjoint) problem for multiple frequencies, 
for which we derive flexible preconditioned Krylov subspace solvers 
specifically for shifted systems. The solvers allow the preconditioner 
to change at each iteration. We analyse the convergence of the solver 
and perform an error analysis when an iterative solver is used for in- 
verting the preconditioner matrices. Finally, we apply our algorithm 
to a challenging application taken from oscillatory hydraulic tomog- 
raphy to demonstrate the computational gains by using the resulting 
method. 



1 Introduction 

Hydraulic tomography (HT) is a method for characterizing the subsurface 
that consists of applying pumping in wells while aquifer pressure (head) re- 
sponses are measured. Using the data collected at various locations, impor- 
tant aquifer parameters (e.g., hydraulic conductivity and specific storage) are 
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estimated. An example of such a technique is transient hydrauhc tomography 
(reviewed in [3] ) . With transient hydrauhc tomography the measured signals 
after a certain distance is typically weak in comparison to ambient noise (i.e., 
disturbances from pumping other other wells, changes in water level in adja- 
cent streams, etc.). Oscillatory hydraulic tomography (OHT) is an emerging 
technology for aquifer characterization that involves a tomographic analysis 
of oscillatory signals. Here we consider that a sinusoidal signal of known 
frequency is imposed at an injection point and the resulting change in pres- 
sure is measured at receiver wells. Consequently, these measurements are 
processed using a nonlinear inversion algorithm to recover estimates for the 
desired aquifer parameters. Oscillatory hydraulic tomography has notable 
advantages over transient hydraulic tomography; namely, a weak signal can 
be distinguished from the ambient noise and by using signals of different 
frequencies, we are able to extract additional information without having to 
drill additional wells. 

Using multiple frequencies for OHT has the potential to improve the qual- 
ity of the image. However, it involves considerable computational burden. 
Solving the inverse problem, i.e. reconstructing the hydraulic conductivity 
field from pressure measurements, requires several application of the forward 
(and adjoint) problem for multiple frequencies. As we shall show in section |5| 
solving the forward (and adjoint) problem involves the solution of shifted sys- 
tems for multiple frequencies. For finely discretized grids, the cost of solving 
the system of equations corresponding to each frequency can be high to the 
extent that it might prove to be computationally prohibitive when many 
frequencies are used, for example, on the order of 200. The objective is to 
develop an approach in which the cost of solving the forward (and adjoint) 
problem for multiple frequencies is not significantly higher than the cost of 
solving the system of equations for a single frequency - in other words, the 
cost should depend only weakly on the number of frequencies. 

Direct methods, such as sparse LU, Cholesky or LDL^ factorization, are 
suited to linear systems in which the matrix bandwidth is small, so that 
the fill-in is somewhat limited. An additional difficulty that direct methods 
pose is that for solving a sequence of shifted systems, the matrix has to be 
re-factorized for each frequency, resulting in a considerable computational 
cost. By contrast, Krylov subspace methods for shifted systems are particu- 
larly appealing since they exploit the shift-invariant property of Krylov sub- 
spaces [28j to obtain approximate solutions for all frequencies by generating 
a single approximation space that is shift independent. Several algorithms 
have been developed for dealing with shifted systems. Some are based on 
Lanczos recurrences for symmetric systems [T^ ^U\; others use the unsym- 
metric Lanczos [9] , and some others use Arnoldi iteration [6l HDl EHJ El [12] . 
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Shifted systems also occur in several other applications such as control the- 
ory, time dependent partial differential equations, structural dynamics, and 
quantum chromodynamics (see [2S] and references therein). Hence, several 
other communities can benefit from advances in efficient solvers for shifted 
systems. The Krylov subspace method that we propose is closest in spirit 
to [12]. However, as we shall demonstrate, we have extended their solver 
significantly. 

Contributions: Our major contributions can be summarized as follows: 

• We have extended the flexible Arnoldi algorithm discussed in [12] for 
shifted systems of the form {K + ajM)xj = b for j = 1, . . . ,nf that em- 
ploys multiple preconditioners of the form {K + tM). Further, analyze 
the convergence of the solver. 

• When an iterative solver is used to apply the preconditioner, we de- 
rive an error analysis that gives us stopping tolerances for monitoring 
convergence without constructing the full residual. 

• Our motivation for the need for fast solvers for shifted systems comes 
from oscillatory hydraulic tomography. We describe the key steps in- 
volved in inversion for oscillatory hydraulic tomography, and discuss 
how the computation of the Jacobian can be accelerated by the use of 
the aforementioned fast solvers. 

Limitations: The focus of this work has been on the computational 
aspects of oscillatory hydraulic tomography. Although the initial results are 
promising, several issues remain to be resolved for application to realistic 
problems of oscillatory hydraulic tomography. For example, we are inverting 
for the hydraulic conductivity assuming that the storage field is known. In 
practice, the storage is also unknown and needs to be estimated from the data 
as well. Moreover, simulating realistic conditions (higher variance in the log 
conductivity field, and adding measurement noise in a realistic manner) may 
significantly improve the error with the addition of information from different 
frequencies. We will deal with these issues in another paper. 

The paper is organized as follows. In section |2j we discuss the Krylov 
subspace methods for solving shifted linear systems of equations based on 
the Arnoldi iteration using preconditioners that are also shifted systems. In 
section [4], we discuss an error analysis when an iterative method is used to 
invert the preconditioner matrices. In section |5| we discuss the basic consti- 
tutive equations in OHT, which can be expressed as shifted linear system of 
equations and discuss the geostatistical method for solving inverse problems. 
Finally, in section [6] we present some numerical results on systems of shifted 
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systems and then discuss numerical results involving the inverse problem aris- 
ing from OHT. We observe significant speed-ups using our Krylov subspace 
solver. 



2 Krylov subspace methods for shifted sys- 
tems 

The goal is to solve systems of equations of the form 

{K + ajM) Xj =b j = 1, . . . , n/ (1) 

Note that aj, for j = 1, . . . , n/ are (in general) complex shifts. We assume 
that none of these systems are singular. In particular, for our application, 
both K and M are stiffness and mass matrices respectively and are positive 
definite but our algorithm only requires that they are invertible. By using a 
finite volume or lumped mass approach, the mass matrices become diagonal 



but this assumption is not necessary. Later, in sections and 5^, we will 
show how such equations arise in our applications. 

Several efficient methods exist for solving the system of equations ([T]), 
which solve for multiple shifts roughly at the cost of a single system. They 
do so by generating a subspace that is independent of the shift and use 
the shift-invariant property of Krylov subspaces ( for a detailed review, see 
[281 section 14.1] and references therein). However, in practice, the number of 
iterations taken can often be large, especially for large matrices from realistic 
apphcations. 

In order to minimize the number of iterations, Meerbergen [12] proposes a 

def 

left preconditioner of the form Kr = K + tM that is factorized and inverted 
using a direct solver. The application of to a vector is, in general, not 
cheap but the spectrum of {K + tM)~^{K + aM) is often more favorable, 
which results in fast convergence of the Krylov methods in just a few iter- 
ations [19]. This form of preconditioning has its roots in solving large-scale 
generalized eigenvalue problems and is known as Cayley transformation [TT] . 
In [21], the authors provide some analysis for choosing the best value of r 
that optimally preconditions all the systems. However, we observed that 
(also, see [12]) using a single preconditioner for all the systems may not yield 
optimal convergence for all systems. In [12], the authors propose a fiexible 
Arnoldi method for shifted systems that uses different values of r resulting 
in different preconditioners at each iteration. This can potentially reduce the 
number of iterations for all the shifted systems. Before we describe our flex- 



ible algorithm in section 2.2 we will derive the right preconditioned version 
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of Krylov subspace method for shifted systems. This serves two purposes - 
it motivates our algorithm, while clarifiying some of the notation. 



2.1 Right preconditioning for shifted systems 

As mentioned earlier, we will review the right preconditioned version of the 
Krylov subspace algorithm for shifted systems. Following the approach in [HI, 
|28] , we solve the system of equations ([T]) using a shifted right preconditioner 
of the form K^ = K + tM 

{K + a,M)K;'x{a,) = h x{a,) = K'x{a,) (2) 

for j = 1, . . . ,n/. 



Algorithm 1 Arnoldi using Modified Gram-Schmidt [24]: 



Given M, h the right hand side. 
Compute Vi = b/P and /3 =^ ||6||2 

dcf 

Choose T and factorize = K + tM 

Define the m + 1 x m matrix Hrn = {hi,k}i<i<m+i,i<k<m- Set Hr, 
for all /c = 1, . . . , m do 

Compute Zk = K~^Vk 
Wk := Mzk 

for alH = 1, . . . , /c do 

hik ■■= wlvi 

Compute Wk := Wk - hikVk 
end for 

hk+i,k ■■= \\wjh- If hk+i,k = stop 

Vk+l = Wk/hk+l^k 

end for 

Define Vm+i = ^i, . . . , Vm+i] and Zm = [zi, . . . , Zm] 



The algorithm proceeds as follows: first, we run m steps of the Arnoldi 
algorithm on the matrix MK~^. This is summarized in algorithm [l| At the 
end of m steps or the Arnoldi process, we have the following relations 

MZm = Vm+lHm (3) 
{K + TM)Zm= (4) 

where, V^Vm = I- Multiplying the first equation by {aj — r) and adding it 
to the second equation gives us 
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(K + ajM)Zm = Vm+i ^ Q + ((T - r)Hr,^ = Vm+iHm{crj; r) (5) 



='-H'm{'7j;r) 

We seek solutions of the form Xm = ZmUm (with zero as the initial guess), 
since by construction, Zm forms a basis for the Krylov subspace }Cm{MK~^ , b), 
where 

ICmiA, h) = span{6. Ah,..., A^-^h} 

By minimizing the residual norm over all possible vectors in }Cm{MK~^,b), 
we obtain the generalized minimum residual (GMRES) method for shifted 
systems, whereas by imposing the Galerkin condition -L /Cm(M_ft'^^, 6), 
we obtain the full orthogonalized method (FOM) for shifted systems. This 
is summarized in algorithm [2j 

Algorithm 2 FOM/GMRES for Shifted Systems 
1: Given matrices K and M, a right hand side b, a E {cxi, . . . , o"„^.} 

^1 1 .1 , def 

2: 
3: 
4: 
5: 
6: 
7: 
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Choose r, build K^- = K + tM and construct preconditioner. 
for all / = 1, 2, . . . , m do 

Run m steps of algorithm [l] to generate Vm+i,Hm and Zm- 

for all J = 1, . . . , n/ do 
If, not converged 

Construct Hjn{(Jj', t) as defined in equation ([5|. 
FOM: 



GMRES: 

yTi.(^j) = min WPei- Hm{,(Jj;r)yJ\\2 

Construct the approximate solution as Xmi^o-j) = ZmUmicrj] 
end for 
end for 



Thus, there is a distinct advantage in using iterative solvers for shifted 
systems; the expensive step of constructing the basis for the Krylov subspace 
is performed only once and using the shift-invariance property of the Krylov 
subspace, the sub-problem for each shift in [2] can be computed at relatively 
low cost. In [T9j, the spectrum of {K + aM){K + tM)~^ was analyzed and 
it was shown that the preconditioner Kr is well suited only for values of 



6 



frequencies a near r. However, the values of a can be widely spread and a 
single preconditioner Kj. might not be a good choice for preconditioning all 



the systems. In section 6.1, we demonstrate an example in which a single 
preconditioner does not satisfactorily precondition all the systems. In [12], 
the authors propose a flexible approach using a (possibly) different precon- 
ditioner at each iteration. We shall adopt this approach. 



2.2 Flexible preconditioning 

We now describe our flexible Krylov approach for solving shifted systems. 
Following Saad [23] and [12], we use a variant of GMRES which allows a 
change in the preconditioner at each iteration. In algorithm [T| we considered 

a fixed preconditioner of the form Kr '= K + tM for a fixed r. Suppose 
we used a different preconditioner at each iteration of the form K + r^M for 
k = 1, . . . ,Tn, then instead of ^ we have, 

{K + TkM)zk = Vk k = l,...,m (6) 

The algorithm is summarized in|4| In this algorithm, in addition to saving Vm, 
we also save the matrix Zm- If at every step in the flexible Arnoldi algorithm 
we use the same value of r, we are in the same position as in algorithm [1} 
We have = [zi, . . . , Zm], H = {hik}i<i<m+i,i<k<m and Kn = [^i, ■■■,Vm] 
which satisfies V^Vm = Im- In addition, we also have the following relations 



KZ„ 



MZrn 



(7) 



where, = diag{ri, 
we obtain for j = 1, . 

{K + ajM)Zm = Vr, 



Tm]- Multiplying ^ by - and adding (|8 



m+l 



(9) 



: H{(Tj]Tm) 



We are now in a position to derive a FOM/GMRES algorithm for shifted 
systems with flexible preconditioning. We search for solutions which are ap- 
proximations of the form Xm(crj) = ZraUmidj)-, which spans the columns of 
Zm- Strictly speaking, spanjZm} is no longer a Krylov subspace. By min- 
imizing the residual norm over all possible vectors in spanjZm}, we obtain 
the flexible generalized minimum residual (FGMRES) method for shifted sys- 
tems, whereas by imposing the Petrov-Galerkin condition -L span{ym,}, 
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l<i<m+l,l<k<m- 



Algorithm 3 Flexible Arnoldi using Modified Gram-Schmidt: 
1: Given M and b the right hand side, Tj,j = l,...,m, vi = b/P and 

P=\\bh 

2: Define the m + 1 x m matrix Hm = {hi^k} 
3: for all = 1, . . . , m do 
4: Solve {K + TkM)zk = Vk 
5: Wk ■■= Mzk 
6: for alH = 1, . . . , do 
7: hik := wlVi 

8: Compute Wk := Wk — hikVk 
9: end for 

10: hk+i,k ■= Ikfclb- If hk+i,k = stop 

11: Vk+1 = Wj/hk+i,k 

12: end for 

13: Define Zm = [zi, . . . , Zm] and Vm+i = ^i, ■ ■ ■ , Vm+i] 



we obtain the fiexible full orthogonalized method (FFOM) for shifted sys- 
tems. This is summarized in algorithm |4| The residuals can be computed 
as 

rmi(^j) = b-{K + ajM)xUaj) (10) 

= Vm+l (/3ei - HmicTk; Tm)ymicrj)) 



2.2.1 Selecting values of Tk 

In equation (|3|, at each iteration we solve a system of the form {K+TkM)zk = 
Vk for k = 1, . . . ,m. This cost can be high if the dimension of the Arnoldi 
subspace m is large and a different preconditioner K + r^M is used at every 
iteration. In practice, it is not necessary to form and factorize m systems 
corresponding to different Tk- In applications, we only need choose a few 
different r^. that cover the entire range of the parameters cTj. The system of 
equations ^ is solved using a direct solver and since only a few values of Tk 
are chosen, the systems can be formed and factorized. Thus, the computa- 
tional cost will not be affected greatly even if the number of frequencies Uf 
is large. 

Let f = {fi, . . . ,fnj,} be the set of values that Tk can take. In other 
words, we take Hp distinct preconditioners. Then, the first mi values of 
are assigned fi, the next m2 values of Tk are assigned f2 and so on. We also 
have m = Ylk=i 
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Algorithm 4 Flexible FOM/GMRES for Shifted Systems 
1: Given matrices K and M, a right hand side b, a E {cxi, . . . , (T„^.} 
2: Choose Tm = diagjri, . . . , r^} and factorize K + t^M. 
3: for all m = 1, 2, . . . do 

4: Run m steps of algorithm [s] to generate Vm+i, Hm+i and Zm- 
5: for all /c = 1, . . . , n/ do 
6: If, not converged 

7: Construct Hm{crk]T.m) as defined in equation ((|9])). 
8: FOM: 

9: GMRES: 

= min ||/3ei - if™(afc;T^)2/m||2 

10: Construct the approximate solution as Xm{,cyk) = ZmUm 
11: end for 
12: end for 



2.2.2 Restarting 

As the dimension of the subspace m increases, the computational and mem- 
ory costs increase significantly. A well known solution to this problem is 
restarting. The old basis is discarded and the Arnoldi algorithm is restarted 
on a new residual. However, for shifted systems, in order to preserve the 
shift-invariant property, one needs to ensure coUinearity of the residuals of 
the shifted systems. For FOM the residuals are naturally coUinear and the 
Arnoldi algorithm can be restarted by scaling each residual by some scalar 
that depends on the shift [12] • For GMRES, the approach used by [TU] 
was extended to shifted systems with multiple preconditioners by [T^. The 
reader is referred to [12] for further details. 

3 Generalized eigenvalue problem and error 
estimates 

We start by computing the approximate eigenvalues and eigenvectors for 
the matrix {K + aM)M'^. Using estimates for approximate eigenvalues 
and eigenvectors we derive expressions for the convergence of the flexible 
algorithms. For convenience, we drop the subscript on the shifted frequency. 
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i.e., use a instead of aj where, j = 1, . . . , -n./. 

Proposition 1. Let Zm, and Vm+i be computed according to algorithm^ 
The eigenpairs of the generalized eigenvalue problem 

Hm{cr; Tm)z = 9HmZ (11) 

then, 9,u = Vm+iHmZ satisfy the Petrov-Galerkin condition J2^ section 
4.3.3] 

{K + aM)M'^u-eu ± span{V^} u e span {Vm+iH^} (12) 

Proof. We first begin by manipulating equations ([T]) and (|8]). Eliminating 
Zm from those equations and adding aVm+iHm to both sides, we have 

Vm+iHm{(rf, Tm) = {K + aM)M-'Vm+iHm (13) 

Now, consider the residual of the eigenvalue calculation for the kth eigenpair, 
where k = 1, . . . ,m is 

r^'^(a) = {K + aM)M-'Vm+lHmZk - OkYm+lHrnZk (14) 

= Ym+lHm{cT]Tm)Zk — OkV„i+lHmZk 

= Vm (Hmicr; Tm)Zk — OkHmZk) — km+l ,mV m+l{Tm + Ok — (y)Kn^k 

= — hm+l,mVm+l {^m + Ok — CT )e*mZk 

From which we can claim that u G span {Kn+i-^m} and {K + aM)M~^u — 
Ou _L span{V^}. In other words, they satisfy the Petrov-Galerkin ( [l2| ) and 
are an approximate eigenpair of {K + aM)M~^ . □ 

Furthermore, we define pk '= \\{K + aM)M~^Uk — OkUk\\2, which is the 
residual of the kth eigenvalue calculations. Small pk implies that the eigen- 
value calculations have converged. It is readily verified that the eigenvalues 
A of KM~^ (and the generalized eigenvalue problem Kx = XMx) are related 
to the eigenvalues A(cr) of {K + aM)M~^ by the relation A(cr) = X + a. 
The importance of the convergence of eigenvalues on the convergence of the 
Krylov subspace solver using FOM can be established by the following result. 

Proposition 2. The residual using FOM satisfies the following inequality 



\rm{(T)h<Y, 



m 

Pk 

k=l 



cr — r„ 



Ok + Tm- cr 



Ok'Wskl 



where, Sk '= e\Z^{H^(3ei) and pk is the residual of the eigenvalue calcula- 
tion, defined above. 



10 



Proof. We start by writing the residual in equation (10) 

Since, for flexible FOM for shifted systems Um = Hm{<^',Tjn)~^ f3ei, the first 
term in the above expression is zero and we have, 



r{a) 



Now, using the generalized eigendecomposition in ( 11 ) Hm{cr; Tm) = HmZQZ ^, 
and using the residual of the eigenvalue calculation in ( 14 ) 



(15) 



k=l 
m 



k=l 



Tm + Ok - (J 



where, Sk = e*k^k H^j3ei. The proof follows from the properties of vector 
norms. □ 

This inequality provides insight into the importance of the accuracy of ap- 
proximate eigenpairs for the convergence of flexible FOM for shifted systems. 
We follow the arguments in [19]. In particular, the residual is very small if 
a ~ Tm, 16*^"^ I, \sk\ or pk are small. We shall ignore the case that a ~ for 
further analysis, i.e. that the shifted system is exactly the preconditioned sys- 
tem. The eigenvalue residual norm pk being small implies that the Ritz values 
are a good approximation to the eigenvalues of {K + (jM)M~^. This implies 
that all the eigenvalues in this interval have been computed fairly accurately. 
We now discuss when |6'^^| is large. When all the values of are equal to r, 
the approximate eigenvalues 9k of KM~^ are related to approximate eigen- 
values \k of the preconditioned system {K + aM){K + tM)~^ by the Cayley 
transformation Therefore, \0^^\ is large only if |Afc + o"| <^ |Afc + r|. The 

term Sk can be rewritten as Sk = elZ^^ H~^V^Vm(3ei = elZ~^H:^^V^b. It 
is readily verified that elZ~^H:^^V^ is orthonormal to all other approximate 
eigenvectors Vm+iH^Zm and thus, Sk can be interpreted as the component 
of the right hand side b in the direction of the approximate eigenvector. In 
this case, the solution has a small component in the direction of b. 

The analysis for the convergence of flexible FOM for shifted systems can 
be extended to flexible GMRES as well. The following result bounds the 
difference in the residuals obtained from m steps using flexible FOM and 
flexible GMRES. 
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Proposition 3. Let Zm, Hm and Vm+i be computed according to algorithm^ 
Further, from algorithm we define the flexible FOM quantities y{^'^{a) = 
H^{a;T^)-^l5ei, residual ri;^^ {a) = V^+i{Pe,-H^{a;T^)y{na)) and flexi- 
ble GMRES quantities yf^^^^'^ (a) = argminy^c'^W (3ei — Hmicr;Tm)y\\2, residual 
rr""'{a) = Vm+i{f3ei-H^{a;Tm)yr""\(^))- Further, assume that H^{a;T^) 
is invertible. We have the following inequality 



< 1 , ll^rr-./,.^ 112 (16) 



r^™(a)||2 " l + \\r]H-*{a-T^)em\\-2 



where, rj = hm+i,m{cr - r^). 



Proof. We begin by the following observation from equation (15) rl!^"^{a) 
-VelaVt"^ and 



Next, we look at the solution to the GMRES least squares problem which 
can be written as 

This can be rewritten as 



In other words, the solution to the GMRES subproblem is a rank-one per- 
turbation of the FOM subproblem. Using the Sherman-Morrison identity 

^ V ' l + \\r]H^*{a;Tm)e^\\i 

Then, the residual difference between FOM and GMRES can be bounded as 



m \ J m 



'l + \\vH-*{a;Tr, 



The inequality ( 16 ) follows from the following observations \\Hm{o", Tm)Hm{(7\ T^) ^ II2 
1 + \\vH-*ia;T^)e^h and {^^^h < \ve*^y'ri^)\. 

□ 
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If ||?7if~*(cr; Tm)em||2 is large, then the difference between the two resid- 
uals can be large. This happens either when rj is large or Hm{o",Tm) is 
close to singular. In this case, GMRES can stagnate and further progress 
may not occur. We now discuss situations in which breakdown occurs, i.e. 
hm+i,m = 0. If hm+i,m 7^ and Hm is full rank, then it can be shown from 
equation ^ that span{MZm} C span{Vm+i} and from equation ([9]), it fol- 
lows that span {{K + (jjM)Zm} C span {Kn+i}. Further, hm+i,m = if and 
only if Xm^cTj) is the exact solution and Hmi^crfTm) is non-singular. The 
argument closely follows |223 and will not be repeated here. 



4 Inexact preconditioning 

We observe that to compute vectors 2;^ for A; = 1, . . . , m in equation ([6]), we 
have to invert matrices of the form K + r^M. When the problem sizes are 
large, iterative methods may be necessary to invert such matrices, resulting 
in a variable preconditioning procedure in which a different preconditioning 
operator is applied at each iteration. More precisely, for = 1, . . . , m, 

~Zk^{K + TkM)-\k Pk = Vk~{K + nM)h (17) 

where, pk is the residual that results after the iterative solver has been ter- 
minated. To simplify the discussion, we assume that the termination criteria 
for the iterative solver is such that \\pk\\2 ^ ^ ll'^fclb = ^, for some e. We also 

=1 

follow closely the approach in [27] . The new flexible Arnoldi relationship is 
now, 

{K + ajM)Zra + Pm = Kn+iiJ™ ((T^- ; T„) j = 1, . . . , (18) 

where, Zm = [zi, ■ ■ ■ , Zm] and Pm = \pi,...,Pm\ and H„,{aj;Tm) is defined 
in equation If we seek approximate solutions of the form Xm{cj) = 

ZmUmio^j), then we can derive expressions for the approximate residuals as 
follows, 

rm(o-i) = h - {K + ajM)Z^ym{cTj) 

= b — Vm+lHm{o'j;Tm)ym{(^ j) + PmUm^f^ j) 

= Vm+1 (/3ei - Hm{(^j] Tra)yra{crj)) +Pmym{<yj) 
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The quantity rm{crj) is the true residual and fm{o'j) is the inexact residual 
that can be computed by ignoring the error due to early termination of the 
iterative solver. In practice, computing the true residual (or even the inexact 
residual) can be expensive. However, in order to monitor the convergence of 
the iterative solver, we need bounds on the true residual. Using such bounds, 
we can derive stopping criteria for the flexible Krylov solvers for shifted 
systems with inexact preconditioning. We start by deriving the difference 
between the true residual and the inexact residual. 

m 

WrmicTj) -rmicrj)\\2 = \\Pmymi0-j)\\2 = || ^ 6^ t/^ ((TjOPfc || 2 (19) 

fc=l 

m 

< ^\(^kym{(^j)\\\Pk\\2 
k=l 

m 

< e^\elym{(Tj)\ =e\\ym{(Tj)\\i 

k=l 



Using this result, the norm of the true residual r^^aj) can be bounded using 
the following relation 

lkm(crj)||2< \\rmi(^j) - rmi(^j)\\2 + \\f.micrj)\\2 (20) 

< e\\ym{crj)\\i + \\/3ei - Hm{,crj]Tm)ym{crj)\\2 



Using this bound on the true residual, we can monitor the convergence 
of iterative solver for each aj. 

Further, using a similar argument to j^T, Proposition 4.1], we can also de- 
rive specialized results for the flexible FOM/GMRES for shifted systems. Let 

r'n^^j) = b-{K+ajM)ZmyTi^j) and rS-'--((T,) = b-{K+ajM)Zmyr'''%^j) 
be the true residual, respectively resulting from the flexible FOM/GMRES 
for shifted systems. We have the following error bounds 



IIO^"(^.)ll2< e\Wri^,)h (21) 
\{Vr^^,Hm{^r,Tjy rrn\2 < e\\HU^f,Tj\\2\\y'r{a,)h 
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5 Application to Oscillatory Hydraulic To- 
mography 



In this section, we briefly review the apphcation of Oscillatory Hydraulic 
Tomography and the Geostatistical approach for solving the resulting inverse 
problem. 



5.1 The Forward Problem 

The equations governing ground water flow through an aquifer for a given 
domain Q with boundary dQ = dQc U dQjy, OVLd fl dVL^ = are given by, 



^ 9^1^ - V ■ (ir(x) V/i(x, t)) = g(x, t), ^en (22) 
at 

/i(x,t)=0, ^edVto (23) 

V/i(x,t) ■n = 0, ^edVtN (24) 

where S's(x) [L"-*^] represents the speciflc storage and -ft'(x) [L/T] represents 
the hydraulic conductivity. In the case of one source oscillating at a fixed 
frequency u [radians/T] , g(x, t) is given by 

g(x, t) = QoSi^ — Xs) cos{ut) (25) 

To model periodic simulations, we will assume the source to be a point source 
oscillating at a known frequency u and peak amplitude Qq at the source 
location x^. In the case of multiple sources oscillating at distinct frequencies, 
each source is modeled independently with its corresponding frequency as in 



(25), and then combined to produce the total response of the aquifer. 

Since the solution is linear in time, we assume the solution (after some 
initial time has passed) can be represented as 

/i(x, t) = 3ft(<l>(x) exp{iujt)) (26) 

$(x), known as the phasor, is a function of space only and contains informa- 
tion about the phase and amplitude of the signal. Assuming this solution. 



the equations (27) in the phasor domain are. 



-V-(ir(x)V$(x))+iu;^,(x)<l>(x) = go5(x-x,), ^eQ (27) 

$(x) =0, X G drio 

V$(x) • n = 0, X G dVlN 
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The differential equation (27) along with the boundary conditions are dis- 
cretized using FEniCS flQ[ [T7t [18] by using standard linear finite elements. 
Solving it for several frequencies results in system of shifted equations of the 
form 



J 



{2i 



where, K and M are the stiffness and mass matrices, respectively, that arise 
precisely from the discretization of (27). 



5.2 The Geostatistical Approach 

The Geostatistical approach (described in the following papers [131 [El E]) 
is one of the prevalent approaches to solve stochastic inverse problems. The 
idea is to represent the unknown field as the sum of a few deterministic low- 
order polynomials and a stochastic term that models small-scale variability. 
Inference from the measurements is obtained by invoking the Bayes' theo- 
rem, through the posterior probability density function which is composed of 
two parts - likelihood of the measurements and the prior distribution of the 
parameters. Let s(x) G M^" be the function to be estimated, here the log 
conductivity, is modeled by a Gaussian random field. After discretization, 
the the field can be written as s ~ A/'(X/3, Q). Here X is a matrix of low- 
order polynomials, /3 are a set of drift coefficients to be determined and Q is 
a covariance matrix with entries Qij = K(xj,Xj), and ■) is a generalized 
covariance kernel 0. The measurement equation can be written as, 

y = h{s)+v, vr^^^{0,R) (29) 

where y G represents the noisy measurements and f is a random vector of 
observation error with mean zero and covariance matrix R. The matrices R, 
Q and X are part of a modeling choice and more details to choose them can 
be obtained from the following references [13]. The operator h : — )• M.^^ 
is the known as the parameter-to-observation map or measurement operator, 
with entries hi that takes the form, 

hi= [ ^{e'^'^5{x-Xi)}dn (30) 

where Xj, is the location of the measurement sensor. 

Following the geostatistical method for quasi-linear inversion [13], we 
compute s and /3 corresponding to the maximum-a-posteriori probability 
which is equivalent to computing the solution to a weighted nonlinear least 
squares problem. To solve the optimization problem, the Gauss- Newton algo- 
rithm is used. Starting with an initial estimate for the field Sq, the procedure 
is described in algorithm [5] 
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Algorithm 5 Quasi-linear Geostatistical Approach 



1: Compute the Ny x Ng Jacobian J as, 



dh 
ds 



(31) 



2: Solve the system of equations, 



/ JkQJ^ + R JkX\ f 
V {JkXf V A 




) 



( 



y - h{sk) + JkSk 




(32) 



3: The update Sk+i is computed by. 



Sk+l — 



(33) 



4: Repeat steps 1 — 3 until the desired tolerance has been reached. (If 
necessary, add a line search). 



Algorithm |5j requires, at each iteration, computation of the matrices QJ]: 
and JkQJk ■ Since the prior covariance matrix Q is dense, a straightforward 
computation of QJ^ can be performed in 0{NyN^). However, for fine grids, 
i.e., when the number of unknowns Ng is large, storing Q can be expensive 
in terms of memory and computing QJ^ can be computationally expensive. 
For regular equispaced grids and covariance kernels that are stationary or 
translation invariant, an FFT based method can be used to reduce the storage 
costs of the covariance matrix Q to 0{Ns) and cost of matrix- vector product 
to 0{NslogNs). For irregular grids, the Hierarchical matrix approach can 
be used to reduce the storage costs and cost of approximate matrix-vector 
product to 0{Ns log Ns) for a wide variety of covariance kernels [25] . Thus, in 
either situation, the cost for computing QJ]^ can be done in O^NyNglogNg) 
and the cost of computing JkQJ^ is 0{NsNy log Ns + NgNy). 

5.3 Sensitivity Matrix computation 

Computing the Jacobian matrix Jk at each iteration is often an expensive 
step. Although explicit analytical expressions for the entries are nearly im- 
possible, several approaches exist. One simple approach is to use finite dif- 
ferences, but this approach is expensive because it requires as many A^ -|- 1 
runs of the forward problem, i.e. one more than the number of parameters 
to be estimated. For large problems, on finely discretized grids the number 
of unknowns can be quite large, so this procedure is not feasible. 
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To reduce the computational cost associated with calculating the sensi- 
tivity matrix we use the adjoint state method (see for example, |29]). This 
approach is exact and is computationally advantageous when the number 
of measurements is far smaller than the number of unknowns. For a com- 
plete derivation of the adjoint state equations for the Oscillatory hydraulic 



tomography, refer to |2]. For the type of measurements described in (30 ), the 
entries of the sensitivity matrix can calculated by the following expression 



dhj 



lUJ- — $ 



dsj 



ds 



3 J 



dK 
^ + — V$ • 

as.- 



(34) 



where '^i, is the known as the adjoint solution. It satisfies the following 
system of equations for each measurement 



~V ■{KVmi) + iujSs^i= -5(x-Xi), 

0, 

V^,;(x)-n= 0, 



x G f2 

X G (9^/5 

X G OVIn 



(35) 



for i = 1, . . . ,Ny. The procedure for calculating the sensitivity matrix can 
thus be summarized as follows. Since (34) is evaluated for all sj for each 



Algorithm 6 Computing Sensitivity Matrix 



1. For a given field s^, solve the forward problem for $. 

2. For each measurement i, solve the forward problem for \E'j 



3. Compute the integral in (34) to calculate the sensitivity. 



measurement i, the adjoint state method requires only Ny + 1 forward model 
solves to compute the sensitivity matrix. Thus, when the number of measure- 
ments is far fewer than the number of unknowns, the adjoint state method 
provides a much cheaper alternative for computing the entries of the Ja- 
cobian matrix. This is typically the case in Hydraulic tomography, where 
having several measurement locations is infeasible because it requires dig- 
ging new wells. 

Further, we realize that equation (34) takes the same form as equa- 
tion ( [28] ) for multiple frequencies. Thus, we can use the algorithms developed 
in section |2| to solve the system of equations (34) for as many right hand sides 



as measurements. It is possible to devise algorithms for multiple right hand 
sides in the context of shifted systems [20l [5] but we will not adopt this 
approach. 
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6 Numerical Experiments and Results 



We present numerical results for the Krylov subspace solvers and its applica- 
tion to OHT. As mentioned before, we use the FEniCS software [T6l[T7l[l8] to 
discretize the appropriate partial differential equations. We use the Python 
interface to FEniCS, with uBLASSparse as the linear algebra back-end. 
For the direct solvers we use SuperLU [7] package that is interfaced by 
Scipy whereas for the iterative solver we use an algebraic multigrid pack- 
age PyAMG [Ij , with smoothed aggregation along with BiCGSTAB iterative 
solver. In the following sections, for brevity, we only present results for FOM 
solver but we observed similar results for the GMRES method as well. 



6.1 Krylov subspace solver 

In this section, we present some of the results of the algorithms that we have 
described in section [2j We now describe the test problem that we shall use for 
the rest of the section. We consider a 2D aquifer in a rectangular domain with 
Dirichlet boundary conditions on the boundaries. For the log-conductivity 
field log i^r(x), we consider a random field generated using an exponential 
covariance kernel fi;(x, y) = 4exp(— 2||x — y||2/iv). Other parameters used 
for the model problem are summarized in table [TJ We choose 200 frequencies 
evenly spaced between the minimum and maximum frequencies, which results 
in 200 systems each of size 90, 601. 



Definition 


Parameters 


Values 


Aquifer length 


L 


(m 


) 


500 


Specific storage 


loj 




(m-1) 


-11.52 


Mean conductivity 


/^( 


log 


K) (m/s) 


-11.02 


Variance of conductivity 


<y^ 


(loj 


iK) 


1.42 


Frequency range 








r 27r 27rl 
1-600' 3 -I 



Table 1: Parameters Chosen For Test Problem 



First, we motivate the need for multiple preconditioners to solve the 
shifted system of equations (28). We begin by looking at the number of 
iterations taken by restarted GMRES without a preconditioner and using a 
single preconditioner. We choose a preconditioner of the form K + rM, for 
three different values of r. For illustration purposes, we use a direct solver 
to invert the preconditioned systems. These values represent the minimum 
frequency, the average frequency and the maximum frequency in the parame- 
ter range. The number of iterations corresponding to restarted GMRES (30) 
for each of the system is computed and displayed in figure [T] We observe 
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that the number of iterations corresponding to the systems increases as the 
frequency of the system decreases. When we use a preconditioner K + rM, 
the systems with frequencies nearby |r| converge rapidly. However, systems 
with frequencies further from |r| converge more slowly, the further away they 
are from the frequency of the preconditioned system |r|. This is consistent 
with the analysis in section [3] and in particular, proposition [2j Thus, no 
single preconditioner effectively preconditions all the systems in the given 
frequency range although, not surprisingly, choosing |t| in the center of the 
frequency range seems to be the best choice. Thus, in order to make the 
iterative method competitive, we consider using multiple preconditioners to 
solve the shifted system (28). 

We choose preconditioners of the form K + r^M, k = 1, . . . ,m, where 
m is the maximum dimension of the Arnoldi iteration. From figure [T| it is 
clear that the systems with smaller frequencies converge slower, so we choose 
the values of \Tk\ that are distributed closer to the origin. In particular, we 
choose the values of such that they are evenly spaced on a log scale in 
the domain co G xl- -^^^ example, for Up = 5, the distribution of |r| 
is illustrated in figure [21 Now, let the possible values that |r| can take be 
labeled as f = {fi, . . . , f„p} where, Up is the number of distinct number of 
preconditioner frequencies. Then, the first mi values of are assigned fi, the 
next m2 values of are assigned f2 and so on. We also have m = Yl2=i 
We pick rrik = m/np. If the algorithm has not converged in m iterations, we 
restart using the method in section 2.2.2 if we are using a direct solver as the 
preconditioner. Else, we recycle the same sequence of preconditioners. We 
implemented both algorithms to invert the preconditioner matrices - using 
direct solver and using an iterative solver which is an algebraic multigrid 
preconditioned BiCGSTAB. Using Up = 5 along with the scheme to choose 
the preconditioner frequencies described above, we observed that the number 
of iterations (and hence, matrix-vector products) in both cases were less than 
40. 

Finally, we present the comparison in terms of the run time of our al- 
gorithm compared to solving each system using a direct solver. The results 
are presented in figure |3| In the plots, 'Direct' implies that every system is 
solved individually by a direct solver. 'Flexible' algorithm uses 5 different 
preconditioners (see figure |2] with a direct solver for inverting the precon- 
ditioners, and solves the FOM subproblem, whereas 'Inexact' uses the pre- 
conditioners and inverts the preconditioners using an iterative solver. We 
see that both the 'Flexible' and 'Inexact' algorithms outperform the 'Direct' 
approach even for a small number of frequencies. A relative tolerance of 
lkm(o'j)||2/lko(<7i)||2 < 10""^° was used as stopping criteria for all systems 
j = 1, . . . ,nf and all systems converged within 40 iterations. 
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log Conductivity 










^ ■ 
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20 40 60 80 100 

X 



-13 
-14 



Number of iterations 





Figure 1: (left) The log conductivity field that we use for the test problem 
with 90, 601 grid points, and (right) iteration count for restarted GMRES (30) 
for unpreconditioned case and also with single preconditioners of the form 
K + tM, with |r| G {||))^)^}- Results indicate that for the particular 
choices made, |r| = ^ which is roughly at the center of the frequency range, 
performs best. 
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Figure 2: Choice of frequencies |r| for the preconditioners K+TkM. Here, for 
illustration, we choose Up = 5 and the distinct values of the preconditioned 
frequencies are evenly space on a log scale in the domain u G 



In figures |4] and [5| we plot the residuals and error as a function of the 
frequency of the system. A direct solver was used for figure |4| whereas an 
iterative solver (algebraic multigrid preconditioned BiCGSTAB) was used in 
figure [5] with a stopping critereon e = 10~^^ (refer to section|4j A relative tol- 
erance of ||rm(o'j)||2/||ro(crj)||2 < 10~^^ was used for all systems j = 1, . . . ,n/ 
and all systems converged within 40 iterations. The behavior of residual and 
the error is quite similar but this is to be expected. 

As discussed in section |4| when an inexact preconditioner is used, the 
flexible Arnoldi relation is no longer exact and the true residual Vmicj) and 



the inexact residual Tmicj) are no longer exactly equal. In fact, the error 
between them can be bounded by the relation (19). In figure, we compare 



the difference between the true residual and the inexact residuals with the 
predicted bound £:||?/m(crj)||i. We see [6] that the bound is fairly accurate. 
The stopping tolerances were chosen to be e = 10~^°, 10"^^, 10^^'^. In fact, 
for larger stopping tolerances for the inner solver, the outer solver did not 
converge. 
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irect solver 
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Figure 3: Comparison of time for the different algorithms. 'Direct' imphes 
that every system is solved individually by a direct solver. 'Flexible' al- 
gorithm uses 5 different preconditioners (see figure [2] with a direct solver 
for inverting the preconditioners, and solves the FOM subproblem, whereas 
'Inexact' uses the preconditioners and inverts the preconditioners using an 
iterative solver. We see that both the 'Flexible' and 'Inexact' algorithms 
outperform the 'Direct' approach even for a small number of frequencies. 



Residual vs. frequency Error vs. frequency 




Figure 4: (left) residual vs frequencies for FOM and GMRES and (right) 
error vs frequency for FOM and GMRES. All systems converged within 40 
iterations. A tolerance of ||rm,(c"j)||2/lko(c"j)||2 < 10~^° was used for all sys- 
tems j = 1, . . . , n/. The preconditioner matrices were inverted using a direct 
solver. 
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Figure 5: (left) residual vs frequencies for FOM and GMRES and (right) 
error vs frequency for FOM and GMRES. All systems converged within 40 
iterations. A tolerance of ||rm(o"j)||2/||ro(crj)||2 < lO^^'^ was used as stopping 
criteria for all systems j = l,...,nf. The preconditioner matrices were 
inverted using a iterative solver with e = 10~^^ as an inner stopping criterion. 



, Residual difference compared to prediction vs. frequency 




Figure 6: Difference in the norm between the true and the inexact residu- 
als, when an iterative solver is used as the preconditioner. Three different 
stopping tolerances were considered e = 10~^°, 10~^^, 10~^^. 
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6.2 Application: Tomographic reconstruction 

The objective is now to determine the random conductivity field K{x) from 
discrete measurements of the head h obtained from several pumping tests 
performed with multiple frequencies. Since the conductivity field needs to 
be positive, so that the forward problem is well-posed, we consider a log- 
transformation s = logK. The "true" field is taken to be that in figure [7] 
which is a scaled version of Franke's function [S]. We choose the covariance 
matrix Q to have entries Qij = K(xj,Xj), corresponding to an exponential 
covariance kernel 

'|x-y||2 



/t(x, y) = exp 



AL 



where, L is the length of the domain. We also choose R = a^I and X = 

The measurements are obtained by taking as the true log conductivity 
field, the field in figure [7} Then, the phasor is calculated by solving equa- 
tions (27) with the source location and measurement locations given in [7| 



The measurements are collected for each frequency and two pieces of informa- 
tion are recored, the coefficients of the sine and cosine terms in equation ([26]). 
The inverse problem is then solved using these measurements. The frequency 
range that is chosen is w G [27r/150, 27r/30]. 

The time for computing the Jacobian is listed in figure |8] For the iterative 
solver, we use the fiexible FOM solver using direct solver as preconditioner. 
The relative stopping tolerance we used was 10"^°. Since the cost to solve 
systems with multiple frequencies is nearly the same as the cost to solve a 
single system, the time for building the Jacobian is, more or less, indepen- 
dent of the number of frequencies. However, when a direct solver is used to 
independently solve the systems for multiple frequencies, the cost for con- 
structing the Jacobian scales linearly with the number of frequencies. This 
results in significant reduction in the cost for solving the inverse problem, 
since constructing the Jacobian is the most expensive part of solving the 
inverse problem. The disparity in the computation times for the Jacobian 
between the direct approach and the iterative procedure is exacerbated fur- 
ther, with larger problem sizes resulting from finer discretizations. 

Finally, we compare the error in the the reconstruction with multiple 
frequencies. Table [2] lists the error in the reconstruction. We report two 
errors - the first being the error in the entire domain, the second being 
the error in the area enclosed by the measurement wells. While increasing 
the number of frequencies does not seem to improve the error in the entire 
domain, it improves the error in the region enclosed by the measurement 
wells. This is the primary motivation for using multiple frequencies in the 
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Figure 7: (left) the true log conductivity field that we use for the synthetic 
inverse problem, and (right) location of measurement wells and the source. 
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Figure 8: Comparison of time taken for different components in the Jaco- 
bian. "Forward" refers to solving the forward problem for multiple frequen- 



cies, equation (27). "Adjoint" refers to solving the adjoint field for multiple 



frequency at each measurement location, equation (35). "Derivatives" refer 



to forming the inner product to form the rows of the Jacobian, equation (34). 
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Figure 9: Comparison of inversion results for log conductivity with nj = 1 
and n/ = 20 frequencies. The errors in reconstruction are reported in|2| 

inversion. However, beyond a point, the addition of frequencies does not seem 
to reduce the error. This might be because there is no additional information 
that is obtained from the addition of measurements with these frequencies, 
and, to further improve estimation accuracy, one would need to introduce 
more stimulation and observation points [2j. 



7 Conclusions 

We have presented a flexible Krylov subspace algorithm for shifted systems of 



the form (28) that uses multiple shifted preconditioners of the form K + tM . 



The values of r are chosen in order to improve convergence of the solver for 
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Total error 


Error within box 


1 


0.2604 


0.0622 


5 


0.2387 


0.0554 


10 


0.2263 


0.0509 


20 


0.2170 


0.0461 



Table 2: error due to the reconstruction. We report two errors - the first 
being the error in the entire domain, the second being the error in the 
area enclosed by the measurement wells. 

all the shifted systems. The number of preconditioners chosen varies based 
on the distribution of the shifts. A good rule of thumb is that the systems 
having shift a will converge faster if there a preconditioner with shift r that is 
nearby a. When the size of the linear systems is much larger, direct solvers 
are much more expensive. In such cases, preconditioning would be done 
using iterative solvers. The error analysis in section |4] provides insight into 
monitor approximate residuals without constructing the true residuals. One 
can naturally extend the ideas in this paper to systems with multiple shifts 
and multiple right hand sides using either block or deflation techniques. 

We applied the flexible Krylov solver to an application problem that ben- 
efited significantly from fast solvers for shifted systems. In particular, oscil- 
latory hydraulic tomography is an emerging technique for aquifer charace- 
terization. However, since drilling observation wells to obtain measurements 
is expensive, oscillatory hydraulic tomography aims to generate more in- 
formative measurements by pumping at different frequencies using the same 
pumping locations and measurement wells. In future studies we aim to study 
more realistic conditions for tomography, including a joint inversion for stor- 
age and conductivity. This would be ultimately beneficial to the practioners. 
We envision that fast solvers for shifted systems would be beneficial for rapid 
aquifer characterization. 
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